Exact solution of a Levy walk model for anomalous heat transport 
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The Levy-walk model is known to provide a good description of anomalous heat conduction in 
one-dimensional systems. In this model the heat carriers execute Levy-walks instead of normal 
diffusion as expected in systems where Fourier's law holds. Here we calculate exactly the average 
heat current, the large deviation function of its fluctuations and the temperature profile of the Levy- 
walk model maintained in a steady state by contact with two heat baths (the open geometry). We 
find that the current is nonlocally connected to the temperature gradient. As observed in recent 
simulations of mechanical models, all the cumulants of the current fluctuations have the same 
system-size dependence in the open geometry. The case of the ring geometry is also discussed. 

PACS numbers: 



Introduction. — Fourier's law of heat conduction (in 
one dimension) states that 



J(x,t) 



dT(x, t) 
dx 



(1) 



where T(x, t), J(x, t) are the local temperature and heat 
current density fields and k is the thermal conductiv- 
ity. Based on results obtained from a large number of 
numerical simulations and various analytical approaches 
it is now believed that Fourier's law is not valid in one 
and two dimensional mechanical systems (in particular 
models where momentum is conserved) and heat con- 
duction is anomalous ^[3]. Anomalous behavior in heat 
conduction is not only a theoretical issue but recently 
of experimental relevance in several low-dimensional ma- 
terials [H [S]. Indicators of the anomaly include — (i) 
in steady states the dependence of the heat current J 
on system size L shows the scaling behaviour J ~ 
with a > 0, (ii) the temperature profiles across systems in 
nonequilibrium steady states arc found to be nonlinear, 
even for very small applied temperature differences and, 
(iii) the spreading of heat pulses in anharmonic chains is 
super-diffusive. 

The microscopic basis of Fourier's law lies in the fact 
that the carriers of heat in a system execute random 
walks. The simplest "derivation" of Fourier's law, from 
kinetic theory, is based on this picture and leads to an 
expression of the conductivity of a material in terms of 
the mean free path £, mean velocity v and specific heat 
capacity c of the heat carriers: k — cv£ (in one dimen- 
sions). The breakdown of Fourier's law thus also implies 
a breakdown of the random walk picture and of the dif- 
fusion equation 



dT(x, t) 
dt 



Kd 2 T(x,t) 
c dx 2 



(2) 



which describes time-dependent heat transfer (assuming 
k has no temperature-dependence). A number of recent 



studies indicate that a good description of anomalous 
heat conduction in one dimensional systems is obtained 
by modeling the motion of the heat carriers as Levy ran- 
dom walks instead of simple random walks [SH5] ■ Numer- 
ical studies show that the spreading of heat pulses and the 
form of steady state temperature profiles can be correctly 
modeled by means of the Levy walk. In another study 
of a model with stochastic dynamics, it has been shown, 
starting from a Boltzmann-equation approach, that the 
temperature satisfies a fractional diffusion equation cor- 
responding to a Levy stable process [ID] . 

While some analytical understanding has been 
achieved [3] it is desirable to further develop the Levy 
walk theory for anomalous heat conductivity so that (i) 
one can use it in the way as one uses Eqs. ( lp l for normal 
diffusion and (ii) one has a clearer idea of the range of 
applicability of the model. In this Letter we present sev- 
eral exact results for steady state heat transport in the 
Levy walk model in one dimension. For setting up the 
steady state we follow the idea in [9] of connecting two 
infinite reservoirs at different temperatures to the system 
and consider a version where space and time are taken 
to be continuous. 

Our exact results provide several interesting physical 
perspective on anomalous heat transport. The analytic 
solution of the particle-density profile exhibits the non- 
linear (and singular) form typical of temperature profiles 
in ID systems [H The steady state current has the 
power law dependence on the system-size, characteris- 
tic of anomalous diffusion. In contrast to ([I]) for normal 
heat transport, it is non- locally connected to the temper- 
ature gradient. We derive the exact cumulant generating 
function of current for the open geometry. Our results 
show that all cumulants of the integrated current have 
the same system-size dependence as the average current. 
This is consistent with recent numerical results of heat 
conduction in hard particle systems [11] and therefore 



2 



strongly suggests that the Levy-walk model gives a good 
description of anomalous heat conduction not only for 
the average current but also for current fluctuations. We 
also show that the cumulant generating function is con- 
sistent with the additivity principle expression [THl 120] • 
For the ring geometry, also, the size dependence of the 
cumulants obtained in simulations can be recovered 
by introducing a size dependent cut-off in the distribu- 
tion of times of the Levy-walk model. 

Levy diffusion on the infinite line. — In the simplest 
description we think of energy in the system as being ef- 
fectively mediated by particles performing Levy walks. 
Each particle carries a single quantum of energy. There- 
fore the local energy density and energy current at any 
point are directly proportional to the particle density and 
current respectively. The precise definition of the Levy 
walk model that we consider here is as follows. For a 
single particle each step of the walk consists in choos- 
ing a time of flight r from a given distribution </>(t) and 
then moving it at speed v over a distance x — vt in 
either direction, with equal probability. Let us define 
P(x, t)dx as the probability that the particle is in the in- 
terval (x, x + dx) at time t. Thus P(x, f) includes events 
where the particle is crossing the interval (x,x + dx). If a 
particle starts at the origin at time t = 0, the probability 
P(x,t) satisfies 



for large t 



P(x,t) = ^(t)6(\x\-vt) 



(3) 



1 /■* 

+ o / dT( t>( T ) [P{x - vr,t - t) + P(x + vt, t - r) 

2 Jo 



where ip( T ) — J°° dr' 4>{t') is the probability of choosing 
a time of flight > r. The Fourier-Laplace transform of 
P(x, t) can be calculated from (|3| in terms of the Laplace 
transform of the distribution of flight times 4>(t) and this 
gives the time dependence of all the cumulants (x n ) c of 
the position x at time t. 

If the first moments of the flight times (r™) were finite, 
the motion would be diffusive. One would then get from 
^ for the first cumulants of x at large t 

(x 2 ) c (r 2 ) (x 4 ) c ^ 3 (r 2 ^ _ 6 (r 2 )(r 3 ) + (r*) _ 



v 2 t 



(r) 



v 4 t 



Here we will consider time-of-flight distributions with a 
power law decay at large time 
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1< /3 < 2 



(for example <j>{r) — (3/t /(l+r/t ) l3+1 which is used for 
the data shown in figure 1). For this range of /3 the mean 
flight time (r) — J dr t 0(t) is finite but (r 2 ) = oo. 
Using again the Fourier-Laplace transform of ^ one gets 



2 A v 2 



(3-/3)(2-/?)/3 (r) 

4 A v 4 
(5-/3) (4 -/3)/3 (r) 



t\ 7 = 3-/3, (6) 
f>+ 2 . (7) 



We see that for I < /3 < 2 the motion is superdiffusive 
[12j [13] . The solution of ([3| corresponds to a pulse whose 
central region is a Levy-stable distribution with a scaling 
x ~ t 1 /^ , has ballistic peaks of magnitude t 1 " 13 at x = 
±vt and vanishes outside this [5]. 

From the evolution equation ^ , one can also see that 
a density fluctuation of wave number k relaxes exponen- 
tially with a time constant t* obtained from the solution 
of J °° rfr0(r) cos(fcwr) exp[r/i*] = 1. For </>(r) of the 
form ([5]) one gets for small k 



1/t* ~ AT(-/3)cos(7r(l-/3/2)) {kvf/{i 



(8) 



whereas when 4>(t) has a finite (t 2 ), the regime is diffu- 
sive with t* ~ fc 2 . 

Lew?/ diffusion in a finite system connected to infinite 
reservoirs. — Let us now turn to the more interesting case 
of the open geometry where the system is a finite segment 
between (0, L) connected on its two sides to reservoirs. 
In addition to the probability density P(x, t) let us also 
define the quantity Q(x, t)dxdt as the probability that a 
particle has precisely landed in the interval (x, x + dx) 
during the time interval (t, t + dt). Note that at any 
given time, a particle could either have landed at a point 
x or could be passing over the point. We need to set 
up the correct boundary conditions required to construct 
a nonequilibrium current carrying steady state. To do 
so we identify the region x < with the left reservoir 
and the region x > L with the right reservoir. We set 
Q(x, t) = Qi for points on the left reservoir and Q(x, f) = 
Q r for those on the right. In the steady state we have 
Q(x,t) = Q(x) and P(x,t) = P(x), and they satisfy 

Q(x)- [ dx' ^-4>(\x - x'\/v) Q(x') 
Jo 2w 

= ^(x/v) + ^{(L-x)/v] , (9) 

dx'^(\x-x'\/v) Q(x') 
2v 

+ ^x{x/v) + ^ X [{L-x)/v] . (10) 



P(x) 



(5) where, as in 0, ip(t) — dr (f>(r) and x(t) 



J°° dTip(r). The problem is closely related to the escape 
probability [M] of a Levy walker on the interval (0, L): 
if H(x) is the probability that a Levy walker starting at 
position x will first hit the left reservoir (i.e. the region 
x < 0) before it hits the right reservoir (i.e. x > L) it is 
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FIG. 1: Plot of Q(x) for different system sizes, for the Levy 
walk with j3 — 1.5, Qi = 1.0, Q r — 0.5. The data are obtained 
by solving Eqs. ( 9|10 1 with discretized space. The distribution 
of flight times is^rj = /3/t a /(l + r/i )' 3+1 . The inset shows 
P(x)/Q(x), which converges to x(0) in the limit of L — > oo. 



easy to see that H(x) satisfies 
H(x) = J dx' ±4>(\x - x'\/v) H(x') + ^(x/v) ,(11) 

and from one can see that 

Q{x) = (Qi-Q r )H(x) + Q r . 

By following the same mathematical steps as in [14] to 
study equations such as ([9| or (11), one can show that, 
in the large L limit, the solution Q(x) of ^ (and H{x) 
of ( 11 )) satisfy 



dx'ipQx - x'\/v) Sgn(x - x')Q'{x') = . (12) 



with Q(0) = Qi and Q(L) = Q r (and H{0) = 1 and 
H(l) = for |TT])) with a solution of §V2§, for a <p(r) 
decaying as in (5^ which satisfies 



>(x) = ~B[x(L~x)f/ 



2-1 



(13) 



Integrating this and imposing the boundary conditions 
Q(0) = Qi and Q(L) — Q r , one obtains the constant 
B = (Qi- Q r ) r(/3)/r(/3/2) 2 L x - (i and therefore Q(x). 



One can also notice that in ( 10 1 the r.h.s. is dominated, 



for large L, by the range \x' — x\ <ti L and therefore 

P(x) = X (0)Q(x) - (t)Q(x) . (14) 

In Fig. 0, we compare numerical results obtained by 
solving Eqs. ( |9|10 ) with the exact results of Eqs. ( 13|14 ). 
The profiles are nonlinear and look similar to those ob- 
served for temperature profiles in ID heat conduction 

[HE]. 

We next discuss the current. The steady state current 
J{x) at position x is given by 



J(x) = 7) / Q(x - x') Sgn(x') ip(\x'\/v) , (15) 

2 J-oo 



which can be interpreted as the difference between the 
flow from left to right and from right to left. After a 
partial integration and using the fact that Q(0) = Qi 
and Q{L) = Q r , one gets 



J(x) 



dx'x(\x-x'\/v) Q'{x') . (16) 



We note that dJ /dx — gives Eq. ( 12 1 and so the current 



is independent of expected. Evaluating the current 



at x = and using Eq. ( 13 ), we get for large L 



AvP T{/3) r(l-f) _ 1 

J^(Ql-Qr) 2 ^ _ 1} r(f) L . 



■0. 
(17) 



From Eq. ^ we then get the relation a = 7— 1, between 
the conductivity exponent of anomalous transport and 
the exponent for Levy-walk diffusion. This relation for 
Levy diffusion was noted in [6], numerically observed in 
ID heat conduction models and a derivation based 

on linear response theory has recently been proposed |16j . 

In the large L limit by using Eq. (14) in Eq. (16) we 
obtain 



J 



v 

"2M 



dx'x(\x-x'\/v)P'(x') . (18) 



This is the analogue of Fourier's Law Eq. ([I]) in the case 
of normal heat conduction and shows that the current is 
nonlocally connected to the temperature gradient. 

Current fluctuations in the open system. — In the rest 
of this letter, we discuss current fluctuations. 

The only contribution to the current through the sys- 
tems is from walkers injected (at rate pl) from the left 
reservoir which end up (either after a non stop flight or 
a non direct flight) into the right reservoir or walkers in- 
jected (at rate pn) from the right reservoir which end up 
into the left reservoir. In the case Qi = 1 and Q r = 
let V be the rate at which walkers, which will end up 
into the right reservoir, are injected from the left reser- 
voir. For general Qi and Q r , because the walkers are 
independent and because they have no preferred direc- 
tion, one has p^ — Q(P and p^ — Q r V. We can then 
write the characteristic function of the integrated current 
Q = J Q dtJ(t) in time t, 

Z{\) = (e XQ ) =n t (e XJ{t}dt ) 

= n t [(l - p L dt - p R dt) + e x p L dt + e- x p R dt] 



where /i(A) = p L (e x - l)+p R (e x -l 



(19) 



Hence we get for the cumulant generating function of Q 

MA) = J Q^ x -V + Qr(e- x -i) ^ 



Ql — Qr 
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where J = (Qt — Q r )'P. We can check that //(A = 0) 



gives the current J. Since in Eq. ( 20 ) the system-size oc 



curs only through the average current J, all cumulants of 
Q have the same size dependence as the average current. 
This is consistent with recent numerical simulations of 
the ID hard-point alternate mass gas [UJ. Not surpris- 
ingly we also note that the following fluctuation theorem 
symmetry relation [17] is satisfied: 



/i(A) = fi[— A - (In Qi - InQr)} 



Current fluctuations in ring geometry. — In the ring 
geometry, the system consists of a fixed number N of 
particles which perform independent Levy walks on a ring 
of length L. In the steady state, the density of particles is 
uniform. As the walkers are independent, the cumulants 
of the integrated current Q are related to those of the 
displacement x(t) of a single walker on the infinite line 
(in the steady state) 



(22) 



where p is the density on the ring. If the walkers perform 
on the ring the same Levy walks as on the infinite line, 
the cumulants of x(t) and therefore those of Q grow, as 
in ( 6|7 |, faster than linearly with time (same exponent 
but a different prefactor as, on the ring, the walker is in 
its steady state rather than starting a flight at t = 0). 

On the other hand if one introduces a cut-off time 
Tj, ~ L s in the distribution </>(t) (for example by ar- 
guing that tl should be of the order of t* , the relaxation 
time corresponding to the shortest wave number on the 
ring k = 2tt/L, implying from ([8| that S = (3; one could 
alternatively argue that, as for the open geometry all the 
length of the flights cannnot exceed the system size and 
therefore 5 = 1). With such a cut-off tl, the cumulants of 
Q would grow linearly in time ( 4|23 l, with an amplitude 
which depends on the system size and on the cumulant 
considered 

(Q 2 )c ^ ^(2_^)5_i , (Q 4 )c ^ ^(4-0)5-3 ^3) 

This cut-off time tl gives therefore a possible explana- 
tion for the behavior seen in simulations on the ring of 
hard-point alternate gas in figure 3 of [IT] , where the cu- 
mulants grow linearly in time with different system size 
dependence {{Q 2 ) c /t ~ L" 5 and (Q 4 ) c /i ~ L 5 ). Then 
one gets from (23) f3 ~ 5/3 and d ~ 3/2 leading through 
(17) to a value a — 1/3 for the anomalous Fourier's law 



of the hard-point alternate gas in the open geometry con- 
sistent with most of the simulations done so far [31 [TTl ITS] 
for this system . How to predict the values of (3 and S 
from first principles in Hamiltoninan systems is of course 
an interesting open question. 

Discussion. — In this work we have studied the Levy 
diffusion model of anamalous heat transport. We have 



computed the average current (17 1, the energy profile 
( 13|14 1 and the large deviation function of the integrated 
current ( 19 1 , in the open geometry, i.e. when the system 
is connected at its two ends to reservoirs. One remark- 
able result is that all the cumulants of the integrated 
current have the same anomalous size dependence. We 
have also proposed a simple possibility to explain the size 
dependence of the cumulants for the ring geometry. 
An interesting observation is that the generating func- 



(21) tion of the integrated current in Eq. ( 20 ) matches exactly 



with the formula obtained from the additivity principle 
(AP) [J5] which gives an expression for /z(A) in terms of 
the the conductivity D and equilibrium current fluctua- 
tions a defined respectively as 



D(Q) = liniAQ^o LJ/AQ 

a(Q) = L lim^ 0o (Q 2 ) /t . 
The expression for fi(X) from AP is 



(24) 



Map(A) = - 



with A 



K 
T 

■Qi 



1 dQ m 

Qr V 1 + 2Ko(Q) 



-i 2 



dQ 



D{Q) 



1 



y/1 + 2Ka(Q) 



.(25) 



(this is a parametric expression: as K varies, fj, and A 
vary). From our exact results for /u(A) we find D = Lp 
and a = L /x"(A = 0) = 2DQ . Using these in 
Eq. ( |25[ ) and after explicitly performing the integrals we 
find hap (A) = £t(A). This result is somewhat surprising 
since the additivity principle is expected normally to hold 
for diffusive systems (here D and a have a L-dependence, 
whereas in usual diffusive systems they don't). This re- 
minds us of the recent numerical study on the validity 
of the additivity principle in mass-disordered harmonic 
crystals where, even in the anomalous regime, the addi- 
tivity principle expression is found to be accurate [20] . 
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